#ifndef EULER_C
#define EULER_C

#include "euler.h"

void euler( double a, double b, /*the segment's borders*/
            double y_a, /*border condition for y*/
            double dy_a, /*border condition for dy/dx*/
            double (*p) (double), /* y'' + p(y)y' + q(y)y = 0 */
            double (*q) (double), /*see the comment above*/
            int N ) /*on how many pieces the segment is split*/
{
    /* Prerequisites:
     * y' = z
     * z' = - p(y)z - q(y)y
     * dx = (b - a) / N
     * y_{k+1} = y_{k} + z_{k}*dx
     * z_{k+1} = z_{k} - p(y_{k})*z_{k}*dx - q(y_{k})*y_{k}*dx
     */

    int i = 0;
    double dx = (b - a)/((double) N),
           y_k = y_a,
           z_k = dy_a,
           x = a,
           tmp = 0;

    /* Create a data file
     * for storing the points.*/
    FILE* data = fopen("data.txt", "ab");

    fprintf(data, "%lf %lf\n", y_a, dy_a);

    /* Iterate through the segment
     * and fill the data.*/
    for (i=0; i<N; i++)
    {
        tmp = y_k;
        y_k = tmp + z_k*dx;
        z_k = z_k - p(tmp)*z_k*dx - q(tmp)*tmp*dx;
        x += dx;
        fprintf(data, "%lf %lf\n", y_k, z_k);
    }

    fclose(data);
}

#endif /*EULER_H*/
